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Hubbard-like Hamiltonians are widely used to describe on-site Coulomb interactions in magnetic 
and strongly-correlated solids, but there is much confusion in the literature about the form these 
Hamiltonians should take for shells of p and d orbitals. This paper derives the most general s, p and 
d orbital Hubbard-like Hamiltonians consistent with the relevant symmetries, and presents them 
in ways convenient for practical calculations. We use the full configuration interaction method to 
study p and d orbital dimers and compare results obtained using the correct Hamiltonian and the 
collinear and vector Stoner Hamiltonians. The Stoner Hamiltonians can fail to describe properly 
the nature of the ground state, the time evolution of excited states, and the electronic heat capacity. 
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I. INTRODUCTION 

The starting point for any electronic structure calculation 
is the Hamiltonian that describes the dynamics of the 
electrons. Correlated electron problems are often formu¬ 
lated using Hamiltonians that involve only a small num¬ 
ber of localized atomic-like orbitals and use an effective 
Coulomb interaction between electrons that is assumed 
to be short ranged, and thus confined to within atomic 
sites. However, there is widespread variation in the liter¬ 
ature over the form chosen for these Hamiltonians, and 
they often do not retain the correct symmetry proper¬ 
ties. Here we derive the most general Hamiltonians of 
this type that are rotationally invariant and describe elec¬ 
trons populating s, p and d orbitals. 

Despite often missing terms, multi-orbital Hamiltonians 
with on-site Coulomb interactions have been used suc¬ 
cessfully to describe strong correlation effects in sys¬ 
tems with narrow bands near the Fermi energy. Ap¬ 
plications include studies of metal-insulator transitions^, 
colossal magnetoresistant manganite^, and d band su¬ 
perconductors such as copper oxides 3 and iron pnictided^. 
The Hubbard Hamiltonian has also been used to sim¬ 
ulate graphene^, where both theory and experiment 
have found magnetic ordering in nanoscale structures 7 ^. 
However, the lack of consistency in the forms used for 


the Hamiltonians is problematic. Indeed, according to 
DagottcP, “the discussions in the current literature re¬ 
garding this issue are somewhat confusing”. 

The Hubbard Hamiltonian 9 is both the simplest and best 
known Hamiltonian of the form we are considering, and 
is valid for electrons in s orbitals. Hubbard was not the 
first to arrive at this form; two years earlier, in 1961, 
Anderson wrote down a very similar Hamiltonian 10 ! for 
electrons in s orbitals with on-site Coulomb interactions. 
Furthermore, in the appendix of the same paper, the 
Hamiltonian was extended to two orbitals, introducing 
the on-site exchange interaction, J. Unfortunately, this 
extra term is not rotationally invariant in spin space. 

In 1966, Rot hPl also considered two orbitals and two pa¬ 
rameters, producing a Hamiltonian that is rotationally 
invariant in spin space. In 1969 Caroli et al 1 -1 corrected 
Anderson’s Hamiltonian for two orbitals, again making it 
rotationally invariant in spin space. However, the Hamil¬ 
tonians of Roth and Caroli et al. do not satisfy rota¬ 
tional invariance in orbital space. This was corrected by 
Dworin and Narath in 197Cp“ who produced a Hamilto¬ 
nian that is rotationally invariant in both spin and orbital 
space. They generalized the multi-orbital Hamiltonian to 
include all 5 d orbitals, although they only included two 
parameters: the on-site Hartree integral, U, and the on¬ 
site exchange integral, J. 

The next contribution to the multi-orbital Hamiltonian 


2 


was by Lyon-Caen and Cyrotf 1 ^, who in 1975 considered 
a two-orbital Hamiltonian for the e g d orbitals. They 
introduced Uq — U = 2 J, where Uo = V aa:aa is the self¬ 
interaction term, U = V a p iOL p is the on-site Hartree in¬ 
tegral, and J = V a f 3 ^ a is the on-site exchange integral, 
important for linking these parameters. Here V is the 
Coulomb interaction between electrons, and a and /3 are 
indices for atomic orbitals. Lyon-Caen and Cyrot also 
included the pair excitation, Jc\ x ^.c\ 1 ^c/2qc/2,t? where 


ct and c are fermionic creation and annihilation opera¬ 
tors. The new term moves a pair of electrons from a 
completely occupied orbital, 2, on site /, to a completely 
unoccupied orbital, 1, on site /. Such double excitations 
are common in the literature 15 18 . 

In 1978, Castellani et a/P 9 wrote down the three-orbital 
Hamiltonian for the t^ g d orbital^l i n a very clear and 
concise way: its form is exactly what we find here for p 
orbital symmetry. They too make note of the equation 
Uo — U = 2 J and include the double excitation terms. 
The next major contribution was by Oles and Stollhoff 2 ^, 
who introduced a d orbital Hamiltonian with three inde¬ 
pendent parameters. The third parameter, A J, repre¬ 
sents the difference between the t^g and the e g d orbital 
exchange interactions; they estimated A J to be 0.15J. 
In our discussion of d orbitals we adopt their notation. 
Unfortunately, their Hamiltonian is not rotationally in¬ 
variant in orbital space. 

In 1989, Nolting et al 22 1 referenced the Hamiltonian pro¬ 
posed by Oles and Stollhoff but discarded the double ex¬ 
citation terms and ignored the A J term. Their Hamilto¬ 
nian will be called the vector Stoner Hamiltonian in this 
paper. Since then, most papers have used a Hamiltonian 
similar to that of Nolting et al ., with only two parame¬ 
ters; see for example references [2314251 The vector Stoner 
Hamiltonian is often simplified by replacing the rotation- 
ally invariant moment-squared operator, m 2 = m m, by 
m 2 , yielding the collinear Stoner Hamiltonian. 

The multi-orbital Hamiltonians derived here consist 
of one-particle hopping (inter-site) integrals and two- 
electron on-site Coulomb interactions. The on-site 
Coulomb interactions for s, p and d orbital symmetries 
are presented in a clear tensor notation, as well as in 
terms of physically meaningful, rotationally invariant (in 
both orbital and spin space), operators that include all 
of the terms from the previous papers as well as addi¬ 
tional terms of the same order of magnitude that have 
not been previously considered. The d orbital Hamil¬ 
tonian presented in section II D| corrects the d orbital 
Hamiltonian proposed by Oles and Stollhoff' 21 by restor¬ 
ing full rotational invariance in orbital space. For an 
s orbital Hamiltonian, the on-site Coulomb interactions 
may be described using one independent parameter. Two 
independent parameters are required for p orbital sym¬ 
metry, and three for d orbital symmetry. Furthermore, 
the method presented can be extended to f and g orbital 
symmetry and generalized to atoms with valence orbitals 
of multiple different angular momenta. 

The use of a restricted basis set and the neglect of the 


two-electron inter-site Coulomb interactions means that 
screening effects are missing. Thus, even though the pa¬ 
rameters that define the s, p and d Hamiltonians appear 
as bare on-site Coulomb integrals, screened Coulomb in¬ 
tegrals should be used when modelling real systems; this 
significantly reduces the values of the parameters. 
Although it is possible to find the on-site Co ulom b inte¬ 
grals using tables of Slater-Condon parameters 26 27 or the 
closed forms for the integrals of products of three spheri¬ 
cal harmonics found by Gaunt 2 -^ and RacatPHl, the point 
of this paper is to remove that layer of obscurity and 
present model Hamiltonians written in a succinct form 
with clear physical meaning. The Slater-Condon param¬ 
eters and Gaunt integrals were not used in the deriva¬ 
tions of the on-site Coulomb interaction Hamiltonians 
presented here; however, they did provide an indepen¬ 
dent check to confirm that the tensorial forms derived 
here are correct. The link between the spherical har¬ 
monics Coulomb integrals, the Slater-Condon parameters 
and the Racah parameters is explained clearly in 39 t The 
transformation to cubic harmonics is straightforward; an 
example of which has been included for the p-orbitals in 
Section HTCl 

We note that RudzikaJ^ 1 has derived rotationally invari¬ 
ant Hamiltonians for p and d orbital atoms. However, his 
derivation is different from ours and his d orbital Hamil¬ 
tonian is partially expressed in terms of tensors; the form 
we present in this paper is written in terms of physically 
meaningful operators. 

The p and d symmetry Hamiltonians are described in 
Section [IT] (see Appendix |B| for full derivations) and com¬ 
pared with the vector Stoner Hamiltonian. In Section 
III] we show the ground state results for both our Hamil¬ 
tonian and the vector Stoner Hamiltonian for p and d 
orbital dimers. We find that, although the results are of¬ 
ten similar, there are important qualitative and quanti¬ 
tative differences. In Section hyi we show the differences 
that emerge when evolving a starting state for a p or¬ 
bital dimer under our Hamiltonian and the vector Stoner 
Hamiltonian. We also calculate a thermodynamic prop¬ 
erty, the electronic heat capacity of a p orbital dimer, 
and obtain significantly different results when using our 
Hamiltonian and the vector Stoner Hamiltonian. The 
eigenstates of all Hamiltonians are found using the full 
configuration interaction (FCI) method as implemented 
in the HANDE codP 


II. THEORY 


A. General form of the Hamiltonian 


We start with the general many-body Hamiltonian of 
electrons in a solid, expressed in second-quantized form 
using a basis of localized one-electron orbitals and ig¬ 
noring terms that do not include electronic degrees of 
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freedom^ 


H ~ ^2 y^^IaJB^Ia.gCjB.cT 

IaJp cr 


IJKLoifixi a £ 


(i) 


Here c\ a a and c Ia a are creation and annihilation op¬ 
erators, respectively, for an electron in orbital a on site 
I with spin a. Upper case Roman letters such as /, J, 
K and L refer to atomic sites, the lower case Greek let¬ 
ters a, /?, % and 7 represent localized orbitals, and a and 
C indicate spin (either up or down). The one-electron 
contributions to the Hamiltonian are encapsulated in the 
tia.jp term that includes the electronic kinetic energy and 
electron-nuclear interaction. The electron-electron inter¬ 
action terms are represented by the Coulomb matrix el¬ 
ements V Ia Jp,K X L 7 - 

We now make two approximations: we only retain 
electron-electron Coulomb interactions on-site and we re¬ 
strict the basis to a minimal set of localized angular mo¬ 
mentum orbitals. Thus we have just one s orbital per site 
for the s-symmetry Hamiltonian, three p orbitals per site 
for the p-symmetry Hamiltonian, and five d orbitals per 
site for the d-symmetry Hamiltonian. The Hamiltonian 
can then be written as follows: 




Iotjpc\a,aCjp,a 


IolJP 



^2 ^2 

I OiPxi 


( 2 ) 


where = V I(xIf3 IxIl . To improve readability, from 

this point onwards we drop the site index, /, from the 
Coulomb integrals and creation and annihilation opera¬ 
tors as we shall always be referring to on-site Coulomb 
interactions. 

The on-site Coulomb interaction, V a p lX1 ^ is a rotation- 
ally invariant tensor if the atomic-like orbitals are angu¬ 
lar momentum eigenstates. To demonstrate this we note 
that applying the rotation operator R to an atomic-like 
orbital <j) a of angular momentum £ gives 


%>«) = Xi 73>4«( r )> ( 3 ) 

P 


where (R) is the (2£ + 1) x (2£ + 1) matrix that corre¬ 
sponds to the rotation R in the irreducible representation 
of angular momentum £. By making the change of vari¬ 
able x = Rr and x' = RE, it is straightforward to show 
that 


V r 


aP,X7 


= / drdr' 


,CaM7<x'( r O0xaM7<x'( r ') 


r - r' 


= (4'a(R)f (4 /5 (R))*V r a ^, xV 4 x (R)4 7 (R), 


(4) 


which demonstrates that V a p :X ^ is indeed a rotationally 
invariant tensor. 


B. The one band Hubbard model: s orbital 
symmetry 

The Hubbard Hamiltonian is applicable only to one band 
models where the one orbital per atom has s symmetry. 
This gives a simple form for the on-site Coulomb inter¬ 
action tensor, Uo = V aa , aa ; as there is only one type of 

orbital, there is only one matrix element. We can use this 

result to simplify the Coulomb interaction part of Eq. © 
and express it in terms of the electron number operator or 
the magnetic moment vector operator, defined as follows: 

™ = XXc7c> ( 5 ) 

a C 

C°CC' a a,C'’ ( 6 ) 

«CC' 

_ _ ( -X _y _£ \ 

^CC' WcC , , cr CC / , a CC , b 

where a are the Pauli spin matrices and the sum over a 
has only one term as there is only one spatial orbital for 
the s case. Hence 

^EE = \ u 0 : ( ? ) 

&Px 7 

or, equivalently, 

^EE V *P,x = -^0 : m 2 ( 8 ) 

OcPx 7 

where we have use the normal ordering operator, ::, to 
remove self interactions. The action of this operator is to 
rearrange the creation and annihilation operators such 
that all the creation operators are on the left, without 
adding the anticommutator terms that would be required 
to leave the product of operators unaltered; if the rear¬ 
rangement requires an odd number of flips, the normal 
ordering also introduces a sign change. For example 

: e ■ = X : 

aPcrC, 

= X 4,c4,<77a,C ( 9 ) 

otPcrC, 

If the normal ordering operator is not used in Eqs. 0 
and © then additional one electron terms have to be 
subtracted to remove the self interaction. 

The mean-field form of this Hamiltonian has been in¬ 
cluded in Appendix | A 1| 


C. The multi-orbital model Hamiltonian: p orbital 
symmetry 

Consider the case where the local orbitals, <a, /?, % and 7 
appearing in Eq. 0 are real cubic harmonic p orbitals 
with angular dependence x/r, y/r , and zjr. In this case 



4 


the rotation matrices, d £ a , a ( R), are simply 3x3 Carte¬ 
sian rotation matrices. This makes Vk/ 3, X7 a rotationally 
invariant fourth-rank Cartesian tensor, the general form 
of which id^U, 


Va/3,x 7 — US a ^Spj -f- JSajfi(3x T ^ ^aft^x 7’ (10) 

where U = V a p ia p, J = V a p^ a and J' = V aa ^p, with 
a j- /3. By examination of the integral in Eq. Q, we 
note that the Coulomb tensor has an additional symme¬ 
try when the orbitals are real (as are the cubic harmonic p 
orbitals), namely V a(3 , xl = V xf3 ^ and V a p iXri = V ai , x p, 
and hence J must be equal to J'. Thus we find 


Va(3,X 7 — U5 ax 5(3 7 + J (S a ^Sp x + ^a/3^x7 )• (ii) 


This recovers the well known equation Uo = U + 2J, 
where Uo = and shows that the most gen¬ 

eral cubic p orbital interaction Hamiltonian is define d by 
exactly two independent parameters—i In Appendix B 
we give a full derivation of Eq. © using representation 
theory. 

In passing we observe that symmetric fourth-rank 
isotropic tensors can be found in other areas of physics, 
such as the stiffness tensor in isotropic elasticity theory^ 


Ciklm — ^ik^lm T T $im$kl)•> 


( 12 ) 


where the vector angular momentum operator for p or¬ 
bitals on site / is defined to be 


e 23ai e 36a)Cg,(jCi3' C7 i (17) 

after 


and e^pa is the three-dimensional Levi-Civita symbol. 
An equivalent expression is 


1 
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EEk 


at 


aft,XIP,C iXX ,c 


aftxi crC 



: n 


2 J- T - 2 

• -- J : m 


J : (flap) : 

aft 


(18) 


where the final term corresponds to on-site electron 
hopping 7 , 


where A is the Lame coefficient and [i is the shear mod¬ 
ulus. 

Transforming Eq. © into a basis of complex spherical 
harmonic p orbitals with angular dependence of the form 
Y^ m (#,</>) is straightforward. We write V rnrn n ^ n i rn m = 

Y< a p xl l *rnd*rn"p l m'x l m'"' 1 V a p, x 7 , where l ma is defined by 

Pm h = ^2 l maPa, (13) 

a 

and has the following values 

( 72 ~72 °\ 

[lma] = ^ 0 1 . (14) 

V 72 ~ 72 0 j 


haft — E! 

a 


(19) 


Eq. (16) embodies Hund’s rules for an atom. By mak¬ 
ing the substitution m = 2S, we see that the spin is 
maximized first (prefactor —2 J) and then the angular 
momentum is maximized (prefactor —\J). 


The mean-field form of the above p orbital Hamiltonian 
is given in Appendix |A 2| 


The result of applying this transformation is 

Ymm" — U$rnm' T J$mm'" $m"m' 

+ (_l)mW J( j_ mm//( j_ m/m///> ( 15 ) 

where m, m', m" and m"' are spherical harmonics indices 
going from -1 to +1, U = Vio,io, and J = Vio,oi- 
Using Eq. © as a convenient starting point, it is 
straightforward to rewrite the on-site Coulomb Hamil¬ 
tonian in terms of rotationally invariant operators? 31 


1 

2 


EEk 


At 


aft,x 7 u a,c r E0,£ C 7>£ C X> cr 


aftxi 


1 

2 


(U-J) 


h 2 : —J : m 2 : —J : L 



(16) 


D. The multi-orbital model Hamiltonian: d orbital 
symmetry 


The on-site Coulomb interaction for cubic harmonic d 
orbitals can be expressed as 


Yap,XI ~ 2 ( USax^Pl 


T ( J T 2^*7 J (fiot'yfifix "t - ^aft^ix) 


48A J ^ ^ ^aab^ftbd^xdt^ta 

abdt 


( 20 ) 






5 


where £ is a five-component vector of traceless symmetric 
3x3 transformation matrices defined as follows 


/_L_ 

2^3 


[£lafc] 


0 
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| 

(° 

0 

1 

2 



0 

0 



0 

0 


(0 

0 

0 

[£3a6] = 

0 

0 

1 

2 

1 

(° 

1 

2 

0 

j 

( 0 

1 

2 

0 

[C4ab] = 

2 

0 

0 

1 

(0 

0 

0 

j 

f 1 

0 

[£5afr] 

0 

- 

1 

2 


0 

1 

2\/3 

0 


0 

0 

1 

•s/3 


( 21 ) 


0 0 


We use the convention that index numbers (1,2, 3,4, 5) 
correspond to the d orbitals (3z 2 — r 2 , zx, yz , xy, x 2 — y 2 ). 
See Appendix B 2 for a derivation of Eq. © using rep¬ 
resentation theory. We have chosen U to be the Hartree 
term between the t^g orbitals, J to be the average of the 
exchange integral between the e g and t^g orbitals, and 
A J to be the difference between the exchange integrals 


for e q and t^ q 


That 


U = V ( 


(zx)(yz),(zx)(yz)i 


J = 


2 {^(zx)(yz),(yz)(zx) T ^{3z 2 —r 2 ){x 2 —y 2 ),{x 2 —y 2 ){3z 2 —r 2 )]^ 

and A J = V^ z 2 _ r 2 ^ x 2 _ y 2 ^^ x 2 _ y 2 ^ Sz 2 _ r 2 ^ - 

V(zx)(yz),(yz)(zx)- This choice of parameters is the 
same as was used by Oles and Stollhofl® Eq. (20) 


generalizes their result and restores rotational invariance 
in orbital space. Note that, in the mean field, the on-site 
block of the density matrix in a cubic solid is diagonal; 
this will cause the mean-field version of Eq. 22) to reduce 
to that of Oles and Stollhofi 21 . However the presence 
of interfaces, vacancies or interstitials will break the 
cubic symmetry make the on-site blocks of the density 
matrix non-diagonal. Therefore the use of the complete 
Hamiltonian is recommended in mean-field calculations 
for systems that do not have cubic symmetry. 

Rewriting Eq. (20) in terms of rotationally invariant op¬ 
erators gives 


Rtot = 


U-^J + 5A j) : n 2 : -1 (J - 6AJ) : m 2 : 


+ (J ~ 6AJ) £ : (h aP f : + Aj : Q 2 : 


aft 


( 22 ) 


where Q 2 is the on-site quadrupole operator squared. 
The quadrupole operator for a single electron is defined 
as 

Qnv — 2 + LjyL/^J — -S^t 2 . (23) 

The quadrupole operator for a system of many electrons 
is obtained by summing the operators for each electron. 


Since each term in the sum acts on the coordinates of 
one electron only, the operators L M and L u appearing in 
that term measure the angular momentum of one elec¬ 
tron only, not the whole system. More detail on the Q 2 
operator is included in Appendi x [C| Eq. (22) is similar 
in form to that found in Rudzikas^T the difference being 
that our Hamiltonian is expressed in terms of physically 
meaningful operators. 

The mean-field form of this Hamiltonian is given in Ap¬ 
pendix |A 3| 


E. Comparison with the Stoner Hamiltonian 


The Stoner Hamiltonian for p and d orbital atoms, as 
generally defined in the literature, has the following form 
for the on-site Coulomb interaction: 25 

Vstoner = \ (u - ^j) ; h 2 : -\j : m 2 (24) 


where J, in this many-body context, is the Stoner pa¬ 
rameter, /. However, in the mean field, the Stoner pa¬ 
rameter should not be taken to be eq ual to J due to the 
self-interaction error; see reference 38 . This form clearly 
breaks rotational symmetry in spin space, which can be 
restored by substituting : m 2 : with : m 2 : 


Trh 2 Stoner 2 2^ 


- 2 t - 2 

: n : —- J : m 


(25) 


In this paper we will compare our Hamiltonians, Eqs. (18) 


and (22), with Eq. (25), which we shall call the vector 
or m^Stoner Hamiltonian. By starting with the vector 
Stoner Hamiltonian and working backwards, it is possible 
to find the corresponding tensorial form of the Coulomb 
interaction: 


Ta/3,x 7 — US ax 5f3ry + J5 a7 5/3 X . 


(26) 


This is very similar to the tensorial form for the p-case 
Coulomb interaction, Eq. (11), except that it is missing 


a term, JS a p5 lx . The lack of this term means that the 
vector Stoner Hamiltonian does not respect the symme¬ 
try between pairs ay and /3y evident from the form of 
the integral in Eq. Q), i.e., ± 1^7°"" As 

we shall see later, this omission changes the computed 
results. 


III. GROUND STATE RESULTS 
A. Simulation setup 

Here we investigate the eigenstates of dimers using both 
our Hamiltonians, Eqs. ( p~8]) and ( [22]) , and the vector 
Stoner Hamiltonian, Eq. (|25|). We employ a restricted 
basis with a single set of s, p and d orbitals for the s, p 
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and d Hamiltonians respectively. We perform FCI cal¬ 
culations using the HANDE 32 computer program to find 
exact eigenstates of the model Hamiltonians. The Lanc- 
zos algorithm was used to calculate the lowest 40 wave- 
functions for our Hamiltonian, more than sufficient to 
find the correct ground state wavefunction, while the full 
spectrum was calculated for the vector Stoner Hamilto¬ 
nian. 

The single particle contribution to the Hamiltonian for 
all models (the hopping matrix) is defined for a dimer 
aligned along the z-axis: in this case it is only non-zero 
between orbitals of the same type, whether on the same 
site or on different sites. The sigma hopping integral, 

| t(j |, was set to 1 to define the energy scale. The other 
hopping matrix elements follow the canonical relations 
suggested by Andersen 39 , and we use the values of Paxton 
and Finnish pp a : pp n = 2 : — 1 for p orbitals, and 
dda : ddir : ddS = —6 : 4 : — 1 for d orbitals. 

The directionally averaged magnetic correlation between 
the two sites has been calculated for the ground state. 
If positive it shows that the atoms are ferromagnetically 
correlated (spins parallel on both sites) and if it is neg¬ 
ative it shows that the atoms are antiferromagnetically 
correlated (spins antiparallel between sites). The correla¬ 
tion between components of spin projected onto a direc¬ 
tion described by angles 0 and 0 , in spherical coordinates, 
is defined as follows 

CW = (: rhl^m 2 ^ :), (27) 

where the :: is the normal ordering operator, used to 
remove self-interaction, and 

™0(t) = ( 2 ^) 

a ££, = cr|^/ cos 0 + cr|^, sin 0 cos 0 + cr|^, sin 0 sin <fi. 
The magnetic correlation averaged over solid angle is 

Cavg = ^(: mi.m 2 :)- (29) 

B. Ground state: p orbital dimer 

The ground states of our Hamiltonian and the vector 
Stoner Hamiltonian have been found to be rather similar 
for 2 , 6 and 10 electrons split over the p orbital dimer, 
but qualitatively different for 4 and 8 electrons. The mag¬ 
netic correlation between the two atoms for 4 electrons is 
shown in Fig. [l] The symmetries of the wavefunctions— 
are indicated on the graph by the notation 2S+1 A^ g : the 
db is only used for L z = 0 (i.e. E states) and corresponds 
to the sign change after a reflection in a plane parallel 
to the axis of the dimer; the u/g term refers to ungerade 
(odd) and gerade (even) and corresponds to a reflection 
through the midpoint of the dimer; A is the symbol cor¬ 
responding to the total value for L z (e.g. L z = 1 is n, 
L z = 2 is A, L z = 3 is <f>, L z = 4 is T); and S is the 



2 4 6 8 10 


U/\t\ 



u/\t\ 


FIG. 1: The magnetic correlation between two p-shell 
atoms, each with 2 electrons, for a large range of 
parameters U/\t\ and J/\t\, where t is the sigma bond 
hopping. The different regions of the graph are labeled 
by the symmetry of the ground state. The top graph is 
generated from our Hamiltonian and the bottom from 
the vector Stoner Hamiltonian. The bottom graph has a 
region with symmetry 3 E“ extending a long way up the 
J axis, which is not present in the ground state of our 
Hamiltonian. The bottom graph also includes a region 
with two degenerate states with symmetries 1 A g and 
1 E+; this degeneracy is broken when our Hamiltonian is 

used. 


total spin. The differences between the ground states 
of our Hamiltonian and the vector Stoner Hamiltonian 
shown in Fig. [l] are as follows: the ground-state wave- 
function of the vector Stoner Hamiltonian has a region 
with 3 E“ symmetry extending far up the J axis, which 
is not present for our Hamiltonian; the region with 1 A g 
symmetry is doubly degenerate for our Hamiltonian, as 
total L z = ±2, whereas for the Stoner Hamiltonian it is 
triply degenerate, as it is also degenerate with the state 
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of symmetry x £+ (this state appears at a higher energy 
for our Hamiltonian). 


C. Ground state: d orbital dimer 


The ground states of our Hamiltonian and vector Stoner 
Hamiltonian have been found to be rather similar for 2, 4, 
6, 14, and 18 electrons split over the d orbital dimer when 
A J is small. The simulations for 10 electrons were not 
carried out as they are too computationally expensive. 
Qualitative differences were found for 8 and 12 electrons; 
for an example see Fig. [2j which shows the magnetic cor¬ 
relation between two atoms in the ground state of the 
d-shell dimer with 6 electrons per atom. From Fig. [2] we 
see that the results for our Hamiltonian with AJ/\t\ = 0.0 
(top graph) and with a small value of AJ/|£| = 0.1 (mid¬ 
dle graph) are qualitatively different from those for the 
vector Stoner Hamiltonian (bottom graph). The vector 
Stoner Hamiltonian makes the and 1 E+ states de¬ 
generate, which is not the case for our Hamiltonian. The 
largest differences are between the AJ/|£| = 0.1 graph 
and the vector Stoner Hamiltonian: new regions with 
symmetry 1 E+, 3 & g and 5 E“ appear, and the region with 
symmetry almost disappears. This shows that the 


inclusion of the quadrupole term in Eq. (22) can make a 
qualitative difference to the ground state. 


IV. EXCITED STATE RESULTS 

Differences between our Hamiltonian and the vector 
Stoner Hamiltonian are also observed for excited states. 
Here we present three examples, one demonstrating ex¬ 
plicitly the effect of the pair hopping, one showing more 
general differences in the excited states through a cal¬ 
culation of the electronic heat capacity as a function of 
temperature, and one showing a difference in the spin 
dynamics of a collinear and a non-collinear Hamiltonian. 
The full spectrum of eigenstates was calculated for all of 
these examples. 


A. Pair hopping 

The pair hopping term, J£ Q(Scr(J , = 

J^ a p : (/W) 2 b is found in both our p and d or¬ 
bital Hamiltonians, but is absent from the vector Stoner 
Hamiltonian. To demonstrate the effect of the pair hop¬ 
ping term we initialize the wavefunction in a state with 
two electrons in the x orbital on site 1 of a p-atom dimer. 
Figure [3] shows the evolution of the wavefunction in time 
using our Hamiltonian and the vector Stoner Hamilto¬ 
nian, with U = 5.0|£| and J = 0.7|£|. For our Hamilto¬ 
nian the pair of electrons in orbital x on atom 1 hops to 
the y and z orbitals on atom 1 (mediated by the pair hop¬ 
ping term) as well as between the two atoms (mediated 



2 4 6 8 10 


U/\t\ 



2 4 6 8 10 


U/\t\ 



U/\t\ 


FIG. 2: The magnetic correlation between two d-shell 
atoms, each with 6 electrons, for a large range of 
parameters U/\t\ and J/\t\, where t is the sigma bond 
hopping. The different regions of the graph are labeled 
by the symmetry of the ground state. The top graph is 
generated using our d-shell Hamiltonian with 
AJ/|t| =0.0; the middle graph is generated using our 
Hamiltonian with a small value ofAJ/|£|=0.1; and the 
bottom graph is generated using the vector Stoner 
Hamiltonian. 



















time(ft/|f|) 



time(ft/|f|) 



k B T/\t\ 


FIG. 4: The electronic heat capacity per atom of a 
dimer with four electrons split over two p-shell atoms 
with U = 5.0|£| and J = 0.7|t|. The solid line is 
generated by solving our Hamiltonian exactly and the 
dashed line by solving the vector Stoner Hamiltonian 
exactly. 


B. Heat capacity 

The heat capacity is calculated from 

^rrnhr^ 2(T>> - (Em)2) - (30) 

where 


FIG. 3: The time evolution under our Hamiltonian 
(top) and the vector Stoner Hamiltonian (bottom) of a 
starting state with two electrons in lx, the x orbital on 
atom 1, of a p-atom dimer with U = 5.01£| and 
J = 0.7|£|. We see that the two electrons in lx are able 
to hop into 1 y and lz when the wavefunction is evolved 
using our Hamiltonian but not when it is evolved using 
the vector Stoner Hamiltonian. 


by the single electron hopping term). The result is that 
there is a finite probability of finding the pair of electrons 
in any of the x, y and z orbitals on either atom. However, 
for the vector Stoner Hamiltonian, the two electrons in 
the x orbital on site 1 are unable to hop into the y and 
z orbitals on atom 1; they are only able to hop to the x 
orbital on atom 2. This means tha there is no possibility 
of observing the pair of electrons in the y or z orbitals 
on either atom. 


(E n (T)) = 


E< e”e x p(-fegr) 
E» ex P(~ fejr) 


(31) 


T is the temperature, ks is Boltzmann’s constant, Si is a 
many-electron energy eigenvalue, and N a is the number 
of atoms. The result for both Hamiltonians as a function 
of temperature for four electrons split over two p-shell 
atoms with U = 5.0|t| and J = 0.7|t| is shown in Fig¬ 
ure [4j There is a qualitative difference between the re¬ 
sults for our Hamiltonian (the solid line) and the vector 
Stoner Hamiltonian (dashed line). The energies of the 
eigenstates of our Hamiltonian are more spread out than 
those of the Stoner Hamiltonian. This is due to the inclu¬ 
sion of the pair hopping term, which causes states with 
on-site paired electrons to rise in energy (by rsj 
The unphysical reduction of the heat capacity to zero as 
the temperature tends to infinity is a consequence of the 
use of a restricted basis set. 


C. Spin dynamics 

Here we show how the collinear Stoner Hamiltonian can 
give rise to unphysical spin dynamics. We consider an 
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S z = +1 triplet with both electrons on one of the two 
p orbital atoms, |T) = |T+i). Written as linear com¬ 
binations of two-electron Slater determinants, the three 
states in the triplet are 

|T + 1 )=-^(-|Lrtlyt> + | 2 *t 2 y f» , 

\T-i) =-t= (~\lx i ly 4,) + \2x i 2 y |)), 

\To) = ~ 4| (|la4 1 V t) + |l^t 1 Vi)) 

+ -h(|2a:4.2yt> + |2st2j / 4.>). (32) 

These states are simultaneously eigenstates of the 
collinear Stoner Hamiltonian, the vector Stoner Hamil¬ 
tonian, and our p-case Hamiltonian. They are degener¬ 
ate, with eigenvalue U — J, for both our p-case Hamil¬ 
tonian and for the vector Stoner Hamiltonian. For the 
collinear Stoner Hamiltonian the states \T+i) and |T_i) 
have eigenvalue U — J whereas state |To) has eigenvalue 
U. This means that the degeneracy of this triplet is bro¬ 
ken in the collinear Stoner Hamiltonian. We now rotate 
|T) in spin space so that the spins are aligned with the 
x-axis, i.e. S x = + 1 , 

l* rot > -\ (\T+i) + |T_!)) + T| To ). ( 33 ) 

|T rot ) is still an eigenstate, with eigenvalue U — J, of both 
the vector Stoner Hamiltonian and our Hamiltonian, but 
it is no longer an eigenstate of the collinear Stoner Hamil¬ 
tonian. The wavefunction |4/ rot ) evolves in time as 

|T tot (r)) = e~^r |^ tot (0)), (34) 


m 2 Stoner and our Hamiltonian 


V. CONCLUSION 

We have established the correct form of the multi-orbital 
model Hamiltonian with on-site Coulomb interactions for 
atoms with valence shells of s, p and d orbitals. The 
methodology used may be extended to atoms with f and 
g shells, and to atoms with valence orbitals of several 
different angular momenta. The results presented show 
that there are important differences between our p- and 
d-shell Hamiltonians and the vector Stoner Hamiltonian. 
The vector Stoner Hamiltonian misses both the pair hop¬ 
ping term, which is present in our p and d orbital Hamil¬ 
tonians, and the quadrupole term, which is present in our 
d orbital Hamiltonian. The pair hopping term pushes 
states with pairs of electrons up in energy, whereas the 
magnetism term pulls states with local magnetic mo¬ 
ments down in energy. The pair hopping term has the 
largest effect on the ground state for p orbitals with 2 
and 4 electrons per atom and for d orbitals with 4 and 
6 electrons per atom, for J/\t\ < 2. This is because the 
number of possible determinants with paired electrons 
on site is large for these filling factors and the low ly¬ 
ing states can be separated based on this. At values of 
J/\t\ > 2 the magnetism becomes the dominant effect 
upon the selection of the ground state and the difference 
between the ground state of our Hamiltonian and that 
of the vector Stoner Hamiltonian becomes small. How¬ 
ever, the differences are rather more pronounced for the 
excited states. This is evidenced by the hopping of pairs 
of electrons between orbitals on the same site, which is 
allowed by our Hamiltonian but not by the vector Stoner 
Hamiltonian, and in differences in the electronic heat ca¬ 
pacity as a function of temperature. We also find clear 
evidence that the collinear Stoner Hamiltonian is inap¬ 
propriate for use in describing spin dynamics as it breaks 
rotational symmetry in spin space. We would expect sim¬ 
ilar problems when using collinear time dependent DFT 
simulations to model spin dynamics. 


f m 2 Stoner and our Hamiltonian 

where r is time. If we take the expectation value of the 
magnetic moment on each site using the collinear Stoner 
Hamiltonian, we find 

<* rot (T)|m lx |* rot (r)} = (v|> rot (r)|m 2 x |vr ot (T)) = cos 

(M/ rot (r)|m ly |^ ot (r)) = (V ot {r)\m 2y \^ (r)) = 0, 
(* rot (T)|m lz |* rot (r)} = (^ rot (r)|m 2z |^ rot (r)) = 0. 

(35) 

In contrast, using the vector Stoner Hamiltonian or our 
p-shell Hamiltonian, the expectation value of the mag¬ 
netic moment is independent of time, equal to (1,0,0) 
in vector format. This demonstrates that calculations 
of spin dynamics using the collinear Stoner Hamiltonian 
can give rise to unphysical oscillations of the magnetic 
moments. 
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Appendix A: Mean-field Hamiltonians 


for which the total energy is 


1. The one band Hubbard model: s orbital 
symmetry 

Application of the mean-field approximation 4 ^ to the on¬ 
site Coulomb interaction part of the Hubbard Hamilto¬ 
nian with s orbital symmetry, Eq. ([7]), yields the follow¬ 
ing, 

Vmf =U 0 h)h - £(5* S c ) s c 2 <^ ; (Al) 


([/-£) (fi) 2 - J -(rh) 2 

J 


+ E 2 NM + M 2 ) 

a/3 

~ E fee*) 

^/W V Z 


a/3a C 

+ ^(CaaCf3()(^a(C/3a) )’ 


(A6) 


where the double counting correction has been included. 


for which the total energy is 


=\u 0 ( <n> 2 - £<5*c c )(cjc ff ) ] , (A2) 

where the double counting correction has been included. 
Equivalently, applying the mean-field approximation to 
Eq. ( 181 ) yields the following, 


3. The multi-orbital model Hamiltonian: d-orbital 
symmetry 

Application of the mean-field approximation®] to the 
model Hamiltonian with d-orbital symmetry, Eq. (22), 
yields the following, 

Vmf = {y “ \ J + 2Al7 ) - l<m).m 


Vmf =|^o ^2(n)n - (m).m - £(5*5444^> ( A3 ) + ^ E ^ n «/3^ n /3« + ^ n a/?) n a/3^ 


for which the total energy is 


=lu 0 ( 2(h) 2 - <m> 2 - £<44X5*54 




(A4) 


where again the double counting correction has been in¬ 
cluded. 


-aj £<<£}<£ 

liv 


~ £ [U(c ] aa c K )c ] K c arj + J ( 5*^44 

a/3aC \ 




+ 48AJ E E ^>ast^>j3tu^>X uv ^l vs ^’/3^X 0 " ) 

a/3'~YX <J C s ^ uv 

(A7) 

for which the total energy is 


2. The multi-orbital model Hamiltonian: p orbital 
symmetry 

Application of the mean-field approximation®] to the 
model Hamiltonian with p orbital symmetry, Eq. (18), 
yields the following, 


Cmf = ( u - (h)h - l(rh).rh 


a/3 


E J d*WV + ^a/3> a . 


a/3 


- E ( u + j 4L4c) 

a(3aC \ 


^a^/3a J 5 


(A5) 


Z?' 

H/ 


Coulomb = 1 (u_\j + 2A Jj (h) 2 - ^<m> 2 

+ E - 2 6AJ (<”aj8)(”/5a> + 

a/3 

+ ^a ./£(<£,><&„> 


a/ 3 cr( ' 

+ ^4L4c)<4cV)) 

+ 24AJ E E ^>otst^> /3tu^>X uv ^> 7 

cx f3'YX (J C stuv 

X (44 7 c><4c£>, ( A8 ) 

where the double counting correction has been included. 


S 7 US 
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Appendix B: Representation theory 


7 ; where 


1. p orbital symmetry 

As a precursor to our treatment of the d case in Appendix 
we use representation theor}@^to derive Eq. ( flTj ) for 
cubic harmonic p orbitals. First we define the follow¬ 
ing irreducible objects: “ 0 ” is a scalar, “ 1 ” is a three- 
dimensional vector, “2” is a 3 x 3 symmetric traceless 
second-rank tensor, “3” is a third-rank tensor consisting 
of three “2”s, and “4” is a fourth-rank tensor consisting of 
three “3”s. The Coulomb interaction V a p^s transforms 
as a tensor product of four “l”s, 1 0 1 0 1 0 1 , but we 
know also that it must be isotropic and that from the 
above irreducible objects only the “ 0 ” is isotropic. 
Objects such as 1 0 1 , represented as Aij, may be ex¬ 
panded just as in the addition of angular momentum: 

1®1 = 0®1®2, (Bl) 


50 — ^ ^ Va(3,xi^ a X^/3'Yi 
X7 

52 = ^ Y^Px = ^ ^ ^^/kX7^ Q /3^X7• (B4) 

X7 xP 

This is equivalent to Eq. <HT where sq = U and 52 = J. 


2. d orbital symmetry 

To find the isotropic five-dimensional fourth-rank tensor 
that describes the on-site Coulomb interactions for d or¬ 
bital symmetry, we find it convenient to map from a five¬ 
dimensional fourth-rank tensor to a three-dimensional 
eighth-rank tensor. We do this by replacing the five¬ 
dimensional vector of d orbitals by the three-dimensional 
traceless symmetric B matrix of the d orbitals: 


where “ 0 ” is a scalar, s = JW AijSij, “ 1 ” is a vector, 

Vi = ^2j k Ajk€ijk, and “ 2 ” is a traceless symmetric ten- 

sor, Mij = |(+ Aji) - \5ij{Yjki A kihi)- The reader 
may be more familiar with combinations of angular mo¬ 
mentum. In the language of angular momentum, the 
transformation properties of a tensor product of two ob¬ 
jects of angular momentum with l = 1 are described by 
a 9 x 9 matrix, which is block diagonalisable into a 1 x 1 
matrix (which describes the transformations of an / = 0 
object), a 3 x 3 matrix (which describes the transforma¬ 
tions of an l = 1 object), and a 5 x 5 matrix (which 
describes the transformations of an l = 2 object). In the 
angular momentum representation the matrix elements 
are complex. Here we are using cubic harmonics and the 
matrix elements are real. Similarly 1 0 1 0 1 0 1 may be 
expanded out as: 

1 0 1 0 1 0 1 = (1 0 1 ) 0 (1 0 1 ), 

= (00 10 2 ) 0(00 10 2 ), 

= (0 0 0 ) 0 (0 0 1 ) 0 (0 0 2 ) 0 (1 0 0 ) 0 (1 0 1 ) 
0(102)0(200)0(201)0(202). (B2) 

The only isotropic, “0”, objects arise from 000, 101 and 
2 0 2 . A general isotropic three-dimensional fourth-rank 
tensor, T^, therefore has three independent scalars that 
are found from symmetric and antisymmetric contrac¬ 
tions of Tijki • We are interested in the form of the p or¬ 
bital on-site Coulomb interactions, V a/ 3 ?X7 , which have an 
additional symmetry, such that they remain unchanged 
under exchange of a and % and under exchange of (3 and 
7 . We therefore only require the symmetric scalars that 
arise from the 0 0 0 and 2 0 2 contractions to describe 
Va/3,X7 : 


Va/3,x 7 — 5oT 52 {^aj^Px T 3 a f3^x 7 ) ’ (B3) 

which is symmetric under exchange of a and x or P and 


B = 


§* 2 -! 0/ 2 + - 2 ) 


xy 


xy 


\y 2 -\{x 2 + z 2 ) 


yz 


yz 

I (z 2 + y 2 ) 


(B5) 


We refer to the elements of B as B a y, and use the sub¬ 
scripts a, 6 , c, d, s, £, u and v as the indices for irre- 
ducibles of rank 2 or higher. The transformation between 
the irreducible B a y and the cubic harmonic d orbitals is 


4>aa (r) = N d ^2^ aa bBab^Rd(r)S(a), (B 6 ) 

ab 


where £ is a traceless, symmetric transformation matrix, 


defined in Eq. (21), Ny = | y ^ is a normalisation factor, 
Ry is the radial function for a cubic d orbital, and S is the 
spin function. The orbital indices <a, /3, x and 7 run over 
the five independent d orbitals, whereas the irreducible 
indices a, 6 , c, d, s, £, u and v run over the three Cartesian 
directions. The mapping between the five-dimensional 
fourth-rank tensor and the three-dimensional eighth-rank 
tensor is as follows 


0qq( r l)^7'( r 2)7^( ri )7^'( r2 ) 


V a p, xy = / dridr 2 


l r l - r 2 | 

= J dridr 2 Arj|i? d (ri)/r^| 2 |i? d (r 2 )/r ^| 2 

E V £* JD* D* p ft t ft 

abed Z-^stuv >>aab ab^>(3cd cd^>x s t st^juv^u 


|ri -r 2 | 


= EE £ Qtab^f3cd^x s t^>'Y uv Vab,cd,st,uv •> 


(B7) 


abed stuv 


where this equation defines the isotropic three- 
dimensional eighth-rank tensor for the on-site Coulomb 
integrals, V a b^ c d,st,uv, and we have dropped the complex 
conjugates on the £ matrices as they are real. The B 
matrix is a three-dimensional traceless symmetric “ 2 ” 
and thus contains 5 independent terms. As we show 
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below, it is possible to represent the isotropic three- 
dimensional eighth-rank tensor as a list of quadruple Kro- 
necker deltas, with 5 independent parameters in general 
and 3 independent parameters for the on-site Coulomb 
interactions. The isotropic three-dimensional eighth- 
rank tensor is then converted back into an isotropic five¬ 
dimensional fourth-rank tensor. For an independent ver¬ 
ification of the form of the isotropic three-dimensional 
eighth-rank tensor see Ref. [44j 

The power of representation theory is that one can follow 
an analogous procedure for shells of higher angular mo¬ 
mentum, using a “3” to represent the seven f orbitals, a 
“4” to represent the nine g orbitals, and so on. One 
can also construct interaction Hamiltonians for atoms 
with important valence orbitals of several different an¬ 
gular momenta. 

We now return to the case of a d shell and explain the 
procedure used to find the general isotropic Coulomb 
Hamiltonian in more detail. We start by represent¬ 
ing the isotropic three-dimensional eighth-rank tensor, 
Vab,cd,st,uv, asa 202020 2 . Proceeding as we did for 
the p case, we write 


2 0 2 0 2 0 2 = (2 0 2 ) 0 (2 0 2 ), 

2(g)2 = 0®l®2®3®4, 
2020202 = (0®1®2®3® 4) 

0 ( 0 ® 1 © 2 ® 3® 4). (B 8 ) 


The terms in Eq. (B 8 ) that can generate a rank 0 (a 
scalar) are 00 0, 10 1, 20 2, 303 and 4 0 4. This im¬ 
plies that a general isotropic three-dimensional eighth- 
rank tensor is defined by five independent parameters. 
We know that, for V a / 3 ?X7 , there exists a symmetry be¬ 
tween the pairs ax and /©y. It follows that V ab ^ c d, s t,uv has 
a symmetry between the pairs ( ab)(st ) and ( cd)(uv ). We 
therefore require only the even contractions, 0 0 0 , 2 0 2 
and 4 0 4, reducing the five parameters to three^ this 
is the same number proposed in Chapter 4 of Ref. [35j 
To get the scalars from V a b, c d,st,uv we have to contract 
it. There are six possible contractions of V a b :C d,st,uv > out¬ 
lined in Table [IJ of which only five are independent. Note 
that one does not contract within the indices a 6 , cd, st or 
uv because B ab is traceless, hence ^2 ab d a bB a b = 0 - By 
examining the diagrammatic contractions from Table [I| 
and maintaining the symmetry between the pairs ( ab)(st ) 
and ( cd)(uv ), we conclude that the relevant contractions 

are A 2 ( A 1 + A 3 ^ A 5 

a c ab,cd,st,uv’> \ ab,cd,st,uv ' ab,cd,st,uv J "> ab,cd,st,uv 

and (Kb,cd,st,uv + A 6 ab,cd,st,uv)- This is evident because 


switching two symmetry-related pairs [that is, switching 
( ab ) and (st) or ( cd ) and ( su )] in diagram 1 yields dia¬ 
gram 3. Similarly, switching two symmetry-related pairs 
in diagram 4 yields diagram 6 . We see that there are four 
contractions but we have only three parameters. This 
implies that out of these four contractions there are only 
three linearly independent contractions. 

We can now write down a general interaction Hamilto¬ 
nian V a b,cd,st,uv us a linear combination of three Kro- 


necker delta products, each multiplied by one of the three 
independent coefficients: 


V r 


ab,cd,st,uv 


^0 ^ab,cd,st,uv 


+ c i ( A « 


ab,cd,st,uv 


■ A 3 i 

^ab^cd^st^uv ) 


■ c 2 A 


ab,cd,st,uv * 


(B9) 


We find by transforming this expression back to 

the five-dimensional fourth-rank representation, as dic¬ 
tated by Eq. (B7). For completeness we include all six of 
the transformed contractions: 


££ A ab,cd,st,uv^aab^/3cd^x s t^ l uv 

abed stuv 

^ v ^aab^/Sba^xst^xts ^a/3^x 7 ’ 
abst 

£ £0 W , uv^>aab^>/3cd^>x s t^ ,r y uv 

abed stuv 

^0 ^aab^xba^Pcd^dc ~ 
abed 

££ a ab : cd : st : uv^>Oiab^^cd^>x s ^ 'yuv 
abed stuv 

^0 ^€atab€'yba€pcd€xdc ~ 
abed 

££ a ab,cd,st,uv^>aab^>/3cd^>x s t^'l uv 

abed stuv 

^0 ^^aab^Pbd^dv^ Xvai 
abdv 

££ a ab,cd,st,uv^aab^/3cd^x s t^ ,r Y uv 

abed stuv 

^0 ^£aab£(3bdfi Xdt^ta’) 
abdt 

££ a ab,cd,st,uv^>aab^>/3cd^>x s t^ l uv 
) 

^ ^ ^-^^aab^xbt^Ptd^da * 


abed stuv 


(B10) 


abdt 


We find that the last three of the transformations of 
the contractions, although rotationally invariant in or¬ 
bital space, are not expressible in terms of Kronecker 
deltas of <a, /?, x an d 7 - This is to be expected as 
there exist terms in V a p, X7 thEit are non-zero and have 
more than two orbital indices, e.g. H( 3 z 2 -r 2 )(xy),(yz)(xz) 

Or F(3 z 2 —r 2 )(xy) ,(xy)(x 2 —y 2 ) • 

We can now write down V a p, xl concisely, 

^c©,X7 = us ax s M + (h + Aj) (5 <27 ^Px ^ol(3$xi) 

- 48A J ^ ZaabZpbdZxdtZ'Yta, (Bll) 

abdt 


where we have defined U to be the Hartree term 
between the t 2g orbitals on-site, U = Vj [ Z x)(yz),(zx)(yz)i 
J to be the average of the exchange integral 
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between the e g and t^g orbitals on-site, J = 

2 {^(zx)(yz),(yz)(zx) “I - ^ r {3z 2 —r 2 ){x 2 —y 2 )^x 2 —y 2 ){3z 2 —r 2 )]i 

and A J to be the difference between the 


exchange integrals for e g and £2 




A J 

V(3z 2 —r 2 )(x 2 — y 2 ),(x 2 — y 2 )(3z 2 —r 2 ) V( Z x)(yz),(yz)(zx) • 

We could equally well have chosen to use the sum of 
the fourth and sixth contractions instead of the fifth 
contraction, although with different coefficients. 


Appendix C: Angular momentum and quadrupole 
operators 

1. Angular momentum 

The angular momentum operators are most naturally ex¬ 
pressed in spherical polar coordinates, even when using 
them to operate on cubic harmonics. Their forms are: 


r . , . , d d 

Ln=i ( sin </>— + cot 0 cos </>—- 
06 00 


d 0 

Ly = i [ - cos 4>— + cot 9 sin (j)— 


f -A 

LZ - ^ A,) 


(Cl) 

(C2) 

(C3) 


where we have set Ti to 1 for simplicity. When applying 
the above operators to the on-site cubic harmonic p or¬ 
bitals (we have dropped the site index for clarity) it is 
straightforward to show that 


LyPj — ^ ^ ^ e yjkPk • 


(C4) 


A general one particle operator O may be expressed in 
terms of creation and annihilation operators as follows: 


O = E ^ocaW^)ci a C 0a ,. 
Oi/3cr<j' 


Substituting Eq. (|C4j) into Eq. (C5) yields 


a/3a 


(C5) 


(C6) 


where fi is a Cartesian direction and a and (3 are cubic 
harmonic p orbitals (x, y or z). The d case is slightly 
mor e com plicated, but is greatly simplified by the use of 
Eq. (B6). Applying the angular momentum operators to 


By findi ng the expectation value and substituting into 
Eq. (C5) we obtain 


4* = *E£ ^yrnk^>akj^>/3jrn^aa^'/3a' (C9) 

jkm a/3a 


2. Quadrupole operator 

We define the quadrupole operator for a single electron 


Q yu — 2 (^LyLy ~b • 


(CIO) 


Here we interpret L g and L v in Eq. (CIO) as compo¬ 


nents of the angular momentum of a single electron. 
The quadrupole operator for a system of N electrons 
is then a sum over contributions from each electron: 
Qyv — Y^i =This makes Q^ v a one-electron 
operator, respresented in second quantization as a lin¬ 
ear combination of strings of one creation operator and 
one annihilation operator. Squaring and tracing the one- 
electron tensor operator Q^ yields a two-electron oper¬ 
ator Q = Q yvQvy 

To find the form of this operator, we start by applying the 
single elect ron version of Q^ to the B matrix, making 
use of Eq. (C7). Starting with L^Ly 


the B matrix reveals that 

Ly\Bjfo) — i ^ ^ {^yjs\Bsk) fi tyks\Bj s ) ) •> (C7) 

s 

with 

Ly \ 4 >ol(j) = INd ^ ^ A ^( cr ) ^2 ^otjk {^yjs\B sk ) fi €yks \Bj s )) . 
jks 

(C8) 


LyL u \B jk ) — iLy E (' e vjs \B sk ) T ^vks \Bj s )) •> 

s 

= {Lj ^ ^ f ^vjs ( e yss'\B s 'k ) fi €yks'\B S s')) 

ss' ' 

T ^vks {p-yjs' |-^s's) T £yss' \Bjs')') ^ i 
= (- 1 ) ^Sjf^Bvk) + 5k^\Bj V ) - 25 vtl \Bjk) 

+ E s^yks' fi ^uks^yjs'^} |-^ss') j 5 (Cll) 

ss' ' 

LyLfi\Bj k ) = ( 1) {^3ju\Bfik) fi $ky\Bjn) 25yy\Bjk) 

fi y ^ (tyjs^vks' fi ^yks^ujs') |-^ss') j 5 (C12) 

ss' ' 

where the last equation was found by exchanging fi and 
v in the previous equ ation . Following the same process, 
the final term in Eq. ( |C10 ) becomes 


yiyLy\B jk ) = 6| B jk ), 


(C13) 


This is not a surprising result as the B matrix contains 
linear combinations of spherical harmonics | l,l z ) wi th l = 
2, and L 2 |/,Z 2 ) = l(l + 1)| l,l z ). Combining Eqs. (Cll), 












14 


(C12) and (C13) we find 


Q ij,v | Bj k ) — {fifjij \B v k) + 

H- Sjw 

Y ^ ^ {^vjs^iiks' Y ^uks^fijs') \B SS ,) j. 
cc / / 


(C14) 


3. The quadrupole operator squared and the 
on-site Coulomb integrals 


The quadrupole operator is of interest to us because one 
of the terms in the d-shell on-site Coulomb interaction 
can be represented in terms of Q 2 . This term is: 


Equation (C5) yields the corresponding operator for a 
system of many electrons: 

Q ^ ^ ^ > {^auk^ftk/j. ^a/ik^ftki^j 


after 

^ ^ y ^>amn^>ftjk^ujm^ l 
mnjk 


pkn J C a<j C ftcT 


(C15) 


We now define the normal ordered quadrupole squared 
as 


: Q 2 := E : ; ■ ( Cl6 ) 

Hv 


By substituting Eq. (C15) into Eq. (C16) and dropping 
the site index for clarity, one obtains the following simple 
formula for the quadrupole squared operator: 


E 

aft 7XCTCT' 


- 36 ax 5 0-t + 9 {^ a -t 5 0 x + 5 <*0 5 xi 


) 


^ ^ v ^ast^fttu^X uv ^7 vs jCacrCft 0 ./Cry a /Cxcr' (^ 17 ) 


■48A J 


E E ^ast^>fttu^>X uv ^l vs ^ CX(J ^ / ft<y' ^x& ’ 

aftx'ycra' stuv 

(CIS) 

which is proportional to the final term in Eq. ( C17| ). The 
other terms in Eq. ( C17| ) also exist elsewhere in the d- 
shell on-site interaction Hamiltonian, so introducing an 
explicit Q 2 term is straightforward: 


-48AJ J2 




it 


ast^fttu^>X uv ^'y vsCacrC ft<j' C 'ycr ,C X cr 


aftxicrcr' stuv 


— ^ AJ : Q : +2 A J ^aa^ftcr'^ftcr'^aa 

aftcrcr' 

- 6AJ (cLc}*'C a <r'Cf>* + cLcL'Cpa'Cpft 

aftcrcr' 

= J : Q 2 : +2AJ : n 2 : +3AJ(: ft 2 : + : m 2 :) 

o 

-6 AJ^: (n a/3 ) 2 :. (C19) 

aft 
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A 


1 

ab,cd,st,uv 



A 


2 

ab,cd,st,uv 


A 


3 

ab,cd,st,uv 


A 


4 

ab,cd,st,uv 


A 


5 

ab,cd,st,uv 


A 


6 

ab,cd,st,uv 


$ad$bcfisv$tu H - ^ac^bd^sv^tu ^ad^bc^su^tv ^ac^bd^su^tv 


fiat&bs&cv&du $as$bt$cv$du H - $at$bs$cu$dv $as$bt$cu$dv 


v^bu^cs^dt T ^au^bv^cs^dt 


^at^bd^cv^su + ^ad^bt^cv^su H - ^at^bc^dv^su H - ^ac^bt^dv^su~\~ 

^at^bd^cu^sv T ^ad^bt^cu^sv T ^at^bc^du^sv H - $ac$bt$du$sv~\~ 
$as$bd$cv$tu T ^ad^bs^cv^tu T ^as^bc^dv^tu T $ac$bs$dv$tu~\~ 
fias$bd$cu$tv $ad$bs$cu$tv H - $as$bc$du$tv T ^ac^bs^du^tv 

^av^bd^ct^su ^ad^bv^ct^su ^av^bc^dt^su ^ac^bv^dt^su~\~ 

fiau&bd&ctfisv H - ^ad^bu^ct^sv T ^au^bc^dt^sv fiac&bu&dt&sv~\~ 
^av^bd^cs^tu T ^ad^bv^cs^tu T ^av^bc^ds^tu T ^ac^bv^ds^tu~\~ 
^au^bd^cs^tv T ^ad^bu^cs^tv H“ ^au^bc^ds^tv T ^ac^bu^ds^tv 


^av^bt^cu^ds T ^at^bv^cu^ds T ^au^bt^cv^ds T ^at^bu^cv ^dsT 
fiav&bs&cu&dt H - 3as&bv&cu&dt T 3au&bs&cv&dt T ^as^bu^cv^dt~\~ 
^av^bt^cs^du T ^at^bv^cs^du T $av$bs$ct$du T ^as^bv^ct^du~\~ 
$au$bt$cs$dv T ^at^bu^cs^dv T ^au^bs^ct^dv T ^as^bu^ct^dv 


TABLE I: Table 


to show the node contractions for a d shell. Every contraction is represented both diagrammatically 

and mathematically. 















